Packages

library(spacetime) # for s-t analysis
library(gstat) # for s-t analysis (krigST)
library(mgcv)
library(mgcViz)
library(leaflet)
library(sp)
library(RColorBrewer)
library(lubridate)

load("/Users/mf/Dropbox (University of Southern California)/Courses/PM569/Fall 2019/Lecture12/STObj3.rda")
STObj4 <- STObj3[, "1993-07-01::1993-07-31"] #subset to July 1993
STObj4[1,] # the data (date, julian date, year, month, day, id, temp, timeindex)
##            julian year month day   id  z timeIndex
## 1993-07-01 728111 1993     7   1 3804 82      1278
## 1993-07-02 728112 1993     7   2 3804 84      1279
## 1993-07-03 728113 1993     7   3 3804 88      1280
## 1993-07-04 728114 1993     7   4 3804 90      1281
## 1993-07-05 728115 1993     7   5 3804 92      1282
## 1993-07-06 728116 1993     7   6 3804 91      1283
## 1993-07-07 728117 1993     7   7 3804 92      1284
## 1993-07-08 728118 1993     7   8 3804 94      1285
## 1993-07-09 728119 1993     7   9 3804 96      1286
## 1993-07-10 728120 1993     7  10 3804 93      1287
## 1993-07-11 728121 1993     7  11 3804 94      1288
## 1993-07-12 728122 1993     7  12 3804 85      1289
## 1993-07-13 728123 1993     7  13 3804 92      1290
## 1993-07-14 728124 1993     7  14 3804 90      1291
## 1993-07-15 728125 1993     7  15 3804 80      1292
## 1993-07-16 728126 1993     7  16 3804 85      1293
## 1993-07-17 728127 1993     7  17 3804 87      1294
## 1993-07-18 728128 1993     7  18 3804 90      1295
## 1993-07-19 728129 1993     7  19 3804 87      1296
## 1993-07-20 728130 1993     7  20 3804 86      1297
## 1993-07-21 728131 1993     7  21 3804 82      1298
## 1993-07-22 728132 1993     7  22 3804 81      1299
## 1993-07-23 728133 1993     7  23 3804 82      1300
## 1993-07-24 728134 1993     7  24 3804 87      1301
## 1993-07-25 728135 1993     7  25 3804 92      1302
## 1993-07-26 728136 1993     7  26 3804 94      1303
## 1993-07-27 728137 1993     7  27 3804 92      1304
## 1993-07-28 728138 1993     7  28 3804 95      1305
## 1993-07-29 728139 1993     7  29 3804 86      1306
## 1993-07-30 728140 1993     7  30 3804 76      1307
## 1993-07-31 728141 1993     7  31 3804 85      1308
pm25 <-read.csv("/Users/mf/Dropbox (University of Southern California)/Courses/PM569/Fall 2019/data/pm25_ca_2018.csv")
pm25_sites<-read.csv("/Users/mf/Dropbox (University of Southern California)/Courses/PM569/Fall 2019/data/pm25_sites.csv")

pm<-merge(pm25, pm25_sites, by='site_id')

pm_pal = colorNumeric(c('darkgreen','goldenrod','brown','brown','brown'),
                        domain=pm$pm25)

  leaflet(pm) %>% 
  addProviderTiles('CartoDB.Positron') %>% 
  addCircles(label=~paste0(pm25, ' ug/m3'), color=~pm_pal(pm25),
             opacity=1, fillOpacity=1, radius=500) %>%
  addLegend('bottomleft', pal=pm_pal, values=pm$pm25,
            title='PM2.5 (ug/m3)', opacity=1)

SpaceTime Variograms

  • We use temperature data and construct separable and joint semivariograms
vv <- variogram(object = z ~ lon + lat, 
                data = STObj4,      # July data
                width = 80,         # spatial bin (80 km)
                cutoff = 1000,      # consider pts < 1000 km apart
                tlags = 0.01:6.01)  # 0 days to 6 days

# st separable variogram
sepVgm <- vgmST(stModel = "separable",
                space = vgm(10, "Exp", 400, nugget = 0.1),
                time = vgm(10, "Exp", 1, nugget = 0.1),
                sill = 20)
sepVgm <- fit.StVariogram(vv, sepVgm)

plot(vv,sepVgm,wireframe=TRUE,all=TRUE)

# st joint variogram (metric, see notes)
metricVgm <- vgmST(stModel = "metric",
                   joint = vgm(10, "Exp", 400, nugget = 0.1),
                   sill = 20,
                   stAni = 100)
metricVgm <- fit.StVariogram(vv, metricVgm)

plot(vv,metricVgm,wireframe=TRUE,all=TRUE)

# st joint variogram (sum metric, see notes)
summetricVgm <- vgmST(stModel = "sumMetric",
                space = vgm(10, "Exp", 400, nugget = 0.1),
                time = vgm(10, "Exp", 1, nugget = 0.1),
                joint = vgm(10, "Exp", 400, nugget = 0.1),
                sill = 20,
                stAni = 100)
summetricVgm <- fit.StVariogram(vv, summetricVgm)

plot(vv,summetricVgm,wireframe=TRUE,all=TRUE)

SpaceTime Kriging

spat_pred_grid <- expand.grid(
                      lon = seq(-100, -80, length = 20),
                      lat = seq(32, 46, length = 20)) %>%
            SpatialPoints(proj4string = CRS(proj4string(STObj3)))
gridded(spat_pred_grid) <- TRUE

## ------------------------------------------------------------------------
temp_pred_grid <- as.Date("1993-07-01") + seq(3, 28, length = 6)

## ------------------------------------------------------------------------
DE_pred <- STF(sp = spat_pred_grid,    # spatial part
               time = temp_pred_grid)  # temporal part

STObj5 <- as(STObj4[, -14], "STIDF")         # convert to STIDF
STObj5 <- subset(STObj5, !is.na(STObj5$z))   # remove missing data


pred_kriged <- krigeST(z ~ lon + lat,       # latitude trend
                       data = STObj5,       # data set w/o 14 July
                       newdata = DE_pred,   # prediction grid
                       modelList = metricVgm,  # semivariogram
                       computeVar = TRUE)   # compute variances

color_pal <- rev(colorRampPalette(brewer.pal(11, "Spectral"))(16))

stplot(pred_kriged,
      main = "Predictions (degrees Fahrenheit)",
      layout = c(3, 2),
      col.regions = color_pal)

pred_kriged$se <- sqrt(pred_kriged$var1.var)
stplot(pred_kriged[, , "se"],
      main = "Prediction std. errors (degrees Fahrenheit)",
      layout = c(3, 2),
      col.regions = color_pal)

## ST GAM

  • We show a separable model using smooths for space and time as well as a space-time interaction model
  • Space is usually defined by a thin plate spline and time as a cubic regression spline. These can be kept as separate components or can be multiplied to create an interactive s-t surface.
  • need the mgcv package
pm$date2<-as.Date(pm$date,"%m/%d/%y")
pm$month <- month(pm$date2)
pm$dow<-wday(pm$date2)
pm$jday = as.numeric(format(pm$date2, '%j'))

# separate space and time
gam1 <- gam(pm25~s(longitude,latitude, k=75,bs="tp") + s(jday,k=15,bs="cr") + s(month,k=8,bs="cc"), data=pm, method="ML")
summary(gam1)
## 
## Family: gaussian 
## Link function: identity 
## 
## Formula:
## pm25 ~ s(longitude, latitude, k = 75, bs = "tp") + s(jday, k = 15, 
##     bs = "cr") + s(month, k = 8, bs = "cc")
## 
## Parametric coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 10.66460    0.07318   145.7   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Approximate significance of smooth terms:
##                          edf Ref.df     F p-value    
## s(longitude,latitude) 65.054  71.36 27.34  <2e-16 ***
## s(jday)               13.848  14.00 88.72  <2e-16 ***
## s(month)               5.783   6.00 21.21  <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## R-sq.(adj) =  0.223   Deviance explained = 22.6%
## -ML =  62973  Scale est. = 91.473    n = 17079
plot(gam1)

vgam1<-getViz(gam1)
plot(vgam1)
## Hit <Return> to see next plot:

## Hit <Return> to see next plot:

## Hit <Return> to see next plot:

gam2 <- gam(pm25~t2(longitude,latitude,jday, k=c(50,6),d=c(2,1),bs=c("tp","cr")) + s(month, bs="cc",k=8), data=pm, method="ML")
pen.edf(gam2) 
## t2(longitude,latitude,jday)rr t2(longitude,latitude,jday)nr 
##                    125.058140                     10.097025 
## t2(longitude,latitude,jday)rn                      s(month) 
##                     80.156901                      5.887925
# 3D plot
vis.gam(gam2)

plot(gam2)

gam3 <- gam(pm25~te(longitude,latitude,jday, k=c(10,6),d=c(2,1),bs=c("tp","cr")), data=pm, method="ML")
pen.edf(gam3) 
## te(longitude,latitude,jday)1 te(longitude,latitude,jday)2 
##                     56.93532                     56.93532
# 3D plot
vis.gam(gam3)